#include "euler.h"
#include  <cmath>
#include  <cassert>

double _gamma_moins_un_ = 0.4;

#define __TORO_1

#ifdef __TORO_1
double _rhol_ = 1.0;
double _vl_ = 0.75;
double _pl_ = 1.0;
double _rhor_ = 0.125;
double _vr_ = 0.0;
double _pr_ = 0.1;
#elif defined __TORO_2
double _rhol_ = 1.0;
double _vl_ = -2.0;
double _pl_ = 0.4;
double _rhor_ = 1.0;
double _vr_ = 2.0;
double _pr_ = 0.4;
#elif defined __TORO_3
double _rhol_ = 1.0;
double _vl_ = 0.0;
double _pl_ = 1000.0;
double _rhor_ = 1.0;
double _vr_ = 0.0;
double _pr_ = 0.01;
#elif defined __TORO_4
double _rhol_ = 5.99924;
double _vl_ = 19.5975;
double _pl_ = 460.894;
double _rhor_ = 5.99924;
double _vr_ = -6.19633;
double _pr_ = 46.0950;
#elif defined __TORO_5
double _rhol_ = 1.0;
double _vl_ = -19.5975;
double _pl_ = 1000.0;
double _rhor_ = 1.0;
double _vr_ = -19.5975;
double _pr_ = 0.01;
#endif

/*----------------------------------------------------------*/

Euler::Euler() : ProblemInterface()
{}

/*----------------------------------------------------------*/

void Euler::computeVisu(Vector& w,const Vector& u) const
{
  double vitesse = u[1]/u[0];
  double pression = _gamma_moins_un_* ( u[2] - 0.5*u[1]*vitesse );
  assert(pression>0.0);
  w[0] = u[0];
  w[1] = vitesse;
  w[2] = pression; 
}

/*----------------------------------------------------------*/

void Euler::initialCondition(Vector& u0, double x) const
{
  if(x < 0.5)
  {
    u0[0] = _rhol_;
    u0[1] = _rhol_*_vl_;
    u0[2] = _pl_/_gamma_moins_un_ + 0.5*_rhol_*_vl_*_vl_;
  }
  else
  {
    u0[0] = _rhor_;
    u0[1] = _rhor_*_vr_;
    u0[2] = _pr_/_gamma_moins_un_ + 0.5*_rhor_*_vr_*_vr_;
  }
}

/*----------------------------------------------------------*/

void Euler::flux(Vector& flux, const Vector& u) const
{
  assert(u[0] > 0.0);
  double vitesse = u[1]/u[0];
  double pression = _gamma_moins_un_* ( u[2] - 0.5*u[1]*vitesse );
  assert(pression > 0.0);
  flux[0] = u[1];
  flux[1] = u[1]*vitesse + pression;
  flux[2] = vitesse * ( u[2]+pression );
}

/*----------------------------------------------------------*/

void Euler::boundaryfluxleft(Vector& f, const Vector& u) const
{
  flux(f, u);
}

/*----------------------------------------------------------*/

void Euler::boundaryfluxright(Vector& f, const Vector& u) const
{
  flux(f, u);
}

/*----------------------------------------------------------*/

double Euler::_computeLambdaMax(const Vector& u) const
{
  double vitesse = u[1]/u[0];
  double pression = _gamma_moins_un_* ( u[2] - 0.5*u[1]*vitesse );
  assert(pression>0.0);
  double c = sqrt( ( _gamma_moins_un_+1.0 )*pression/u[0] );

  return std::max( std::max( fabs(vitesse+c), fabs(vitesse-c) ), fabs(vitesse) );
}

/*----------------------------------------------------------*/

void Euler::interiorflux(Vector& f, const Vector& ul, const Vector& ur) const
{
  int n = f.size();
  Vector fl(n), fr(n);

  flux(fl, ul);
  flux(fr, ur);

  // centré
  f.equal(0.5, fl);
  f.add(0.5, fr);

  // diffusion LF
  double ll = _computeLambdaMax(ul);
  double lr = _computeLambdaMax(ur);
  double diff = std::max(ll, lr);

  for(int i = 0; i < n; i++)
  {
    f[i] += 0.5*diff*( ul[i]-ur[i] );
  }
}

